Rim curvature anomaly in thin conical sheets revisited 

Jin W. Wang 

The James Franck Institute and The Department of Physics, 
The University of Chicago, 929 East 57th Street, Chicago, Illinois 60637 

(Dated: January 20, 2013) 

Abstract 

This paper revisits one of the puzzling behaviors in a developable cone (d-cone), the shape obtained by 
pushing a thin sheet into a circular container of radius R by a distance i] The mean curvature was 
reported to vanish at the rim where the d-cone is supported 02j]. We investigate the ratio of the two principal 
curvatures versus sheet thickness h over a wider dynamic range than was used previously, holding R and 
r] fixed. Instead of tending towards 1 as suggested by previous work, the ratio scales as {h/Ry^^. Thus 
the mean curvature does not vanish for very thin sheets as previously claimed. Moreover, we find that the 
normalized rim profile of radial curvature in a d-cone is identical to that in a "c-cone" which is made by 
pushing a regular cone into a circular container. In both c-cones and d-cones, the ratio of the principal 
curvatures at the rim scales as {R/h)^/"^ F / {Y R^), where F is the pushing force and Y is the Young's 
modulus. Scaling arguments and analytical solutions confirm the numerical results. 

PACS numbers: 46.70.De, 68.55.-a, 46.32.+X 
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I. INTRODUCTION 



When a piece of thin sheet like paper crumples, it develops a network of two types of sharp 
structures: straight stretching ridges and pointlike vertices. Since thin sheets and membranes are 
very common in both natural and man-made structures at almost all length scales, crumpling has 



applications in a broad range of systems, such as graphene sheets yD, carbon nanotubes 



virus capsids [6J, pollen grains [jTi], polymerized membranes |18|], and leaves The pointlike 
vertex singularity has mostly been studied via the simple realization known as the developable 
cone or d-cone, shown in Figure [fla). A d-cone is the shape created by pushing the center of a 
thin sheet into a circular container of radius i? by a distance t] yj]. d-cones, stretching ridges, and 
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crumpling in general have been studied extensively 
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2911 ■ In certain cases it is possible to set 



bounds on the energy of singular structures [|30l Bill , but to our knowledge such bounds have not 



been established for d-cones. Furthermore, while the scaling properties of stretching ridges have 
been determined analytically and numerically [[lll[l2 , 191. several phenomena in d-cones are still 
beyond our understanding Villn . One mysterious behavior of d-cones is the reported vanishing of 
mean curvature at the rim region where a d-cone is supported 

The deformation of a d-cone can be characterized by the deflection e = r]/R and its Young's 
modulus is denoted by Y. r and 9 are defined as radial and angular components in the material 
coordinate system. Crr and Cgg are the radial and azimuthal curvature respectively. In principle 
the shape of a d-cone is governed by the Foppl-von Karman equations [|l3L 13211. whose analytical 



solution is known only in a few special cases [|33 



3411 . Unfortunately, a d-cone is not one of these 



cases. However, it is energetically much cheaper for thin sheets to bend than to stretch, so a d- 
cone is asymptotically unstretched, and thence developable, except in the core region where it is 



25|1 



pushed. Assuming that a d-cone is unstretched almost everywhere, Cerda and Mahadevan [|15 , 
described the deformation in terms of the classical Elastic a of Euler BSll and obtained the shape of 
the d-cone by minimizing the bending energy. Numerical study n2m has confirmed their finding. 

Under the assumption that a d-cone is unstretched almost everywhere, the shape of the d-cone 
has zero radial curvature Crr except in the core region, but the real shape must have nonzero Crr 
at the rim to balance the normal force from the edge of the container[|2i |36|l. Liang and Witten 
M\ reported a striking feature that within numerical accuracy, as the thickness of the sheet went 
to zero, the radial curvature Crr and the azimuthal curvature Cee were nearly equal and opposite 
at the rim, so that the mean curvature, defined as {Crr + Coe)/2, nearly vanished there. They 



2 



(a) 



Q 



(b) 




V 



\ 



Figure 1 : (Color online) (a) A typical simulated d-cone formed by pushing the center O of a hexagonal 
elastic sheet (Equations [T] and ^ against a circular container with concentrated force F. The red solid 
line shows the rim of the container. The side length of the sheet £ is 1.77 times the container radius R. 
The thickness h = R/866, and the displacement of point O is a tenth of R. For clarity, the vertical 
scale is expanded by about six times, (b) A variable lattice for simulating d-cones. It shows the material 
coordinates of the lattice points used in the simulation. The local lattice point densities at the rim and in 
the center are about 4.4 and 3.2 times the overall average point density, respectively. The average distance 
between adjacent lattice points is about R/88. 

found that the feature was independent of container radius R, thicknesses of the sheet h, and 
deflection e, but the circular symmetry of the container was indispensable. Since nonzero Crr 
entails stretching, it should be necessary to consider the stretching in the rim region in order to 
understand this vanishing of mean curvature feature DZfl. 

In the limit of thin sheets, the cost of stretching becomes prohibitive and the rim region with 
nonzero stretching should shrink to zero. In this limit, it seems plausible to treat the rim region 
as a boundary layer sandwiched by regions where the Elastica approach can still be applied. This 
type of boundary layer phenomenon appears in a wide variety of systems, such as the Pogorelov 



ring ridge [|32 



3711 . the "minimal ridge" [|l2ll . and more recent work Da 



38 



Judging from its characteristics, this vanishing mean curvature phenomenon is nonlocal and 
purely geometric; thus, researchers have investigated the connection between this phenomenon 
and another nonlocal geometric constraint on surfaces, i.e. the Gauss-Bonnet theoremls^ 401. 



This theorem requires that the sum of the integral of Gaussian curvature within a surface and the 
integral of geodesic curvature along its boundary remain a constant. For a d-cone, the integral 
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of geodesic curvature along the boundary can be assumed a constant, so net Gaussian curvature 
in one region needs to be balanced by another region(s) of opposite Gaussian curvature. Since 
zero mean curvature at the rim means negative Gaussian curvature, researchers have suspected the 
negative Gaussian curvature at the rim is necessary to balance a net positive Gaussian curvatures in 
the core region [13 6ll . However, the integral of (negative) Gaussian curvature near the rim is almost 
completely compensated by that of two adjacent bands. Thus, the Gauss-Bonnet theorem offers 
no obvious explanation of the nonlocality implicit in the vanishing mean curvature phenomenon. 

This paper investigates the ratio of the two principal curvatures | Crr/Coe \ at the rim versus sheet 
thickness h over a wider dynamic range than was used previously in Ref . [|2|] , holding the deflection 
e and the container radius R fixed. The numerical models are specified in Sec. II. In Sec. Ill, we 
describe our numerical findings in detail. Instead of tending towards 1 as required by the vanishing 
of mean curvature, the ratio \Crr/Cee\ at the rim goes below 1 and scales as (h/RY^^. To better 
understand this power law, we study d-cone's close cousin "c-cone" which is made by pushing 
a regular cone into a circular container, as seen in Figure a) . C-cones are simpler structures 
than d-cones. In a c-cone, deflection e, pushing force F, and thickness h are independent degrees 
of freedom. In both c-cones and d-cones, we find that oc {R/hf/'^F/{YR^) for fixed 

deflection e. Moreover, given the same h, R, and e, the normalized rim profiles of radial curvature 
are identical in a c-cone and a d-cone. To put these numerical findings on firmer grounds, in Sec. 
IV we give scaling arguments for both c-cones and d-cones. General solutions for symmetrically 



loaded conical shells are available 141 



4211 . and we use the proper boundary conditions to get the 



analytical solutions for c-cones. Both the scaling arguments and the analytical solutions confirm 
the numerical results. In Sec. V, we discuss the implications of our findings. 



II. NUMERICAL METHODS 



In principle, one could use standard finite element softwares such as Abaqus 1431] to simulate 
thin sheets, but we find them difficult to adapt to the ultra-thin asymptotic behavior of d-cones that 
we want to study. We use an extended Seung-Nelson model 11191 14411 to cope with the singularity 
at the center of a d-cone and to better study the interaction between the sheet and the supporting 
container. The extended model simulates an elastic sheet by a triangular lattice with variable lattice 
spacing, so it has more adaptability than the original Seung-Nelson model ll45n which dictates a 
uniform lattice. This allows the extended model to concentrate lattice points as needed in regions 
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Figure 2: (Color online) (a) A typical simulated c-cone formed by pushing the center of a regular cone into 
a circular container with concentrated force F. The opening angle of the cone is 168.58°, which translates 
to a deflection of 0.10. The thickness of the elastic sheet (Equations [T] and |2]) h = R/866, and the height 
of the cone is about R/6. For clarity, the vertical scale is expanded by about ten times, (b) A schematic 
diagram of the local deformation near the rim along a meridian of the cone. The depth of the deformation 
is d and width is b. For any point Q, its tangent vector along the meridian is denoted by t, and the deviation 
of t from the original unperturbed meridian is defined as cj). 

of strong gradients. Our deformed lattice is shown in Figure [TJb). 

The total elastic energy of the sheet is the sum of stretching and bending energies on each 
triangle. On an arbitrary triangle the strain tensor and curvature tensor are assumed constant. The 
strain tensor 7 is determined from the changes in the edge lengths, as seen in Figure |3l and the 
curvature tensor C from the dihedral angles between the given triangle and its three adjoining 
triangles as shown in Figure |4l The specific transformation formulas are derived in Appendix A. 
Once we know the strain and curvature tensors, we obtain the corresponding energy densities Es 



and Eb via the conventional equations of elasticity [1321 . 



fiY 

= 2(r^t(Tr(7))' + 2{u - l)Det(7)], (1) 
EB = ^«:(Tr(C))2 + «:GDet(C), (2) 

where u is the Poisson's ratio, k = Yh^/ (12(1 — z/^)) is the bending rigidity, and is the Gaussian 
bending rigidity. The total elastic energy of the sheet is taken as the sum of the energy density for 
each triangle times its undeformed area. 



A variable lattice for a d-cone is created in two stages First a uniform triangular lattice 

of spacing a is used to span the desired area. The lattice spacing a is the distance between adjacent 
lattice points. Then this lattice is mapped to the desired nonuniform lattice e.g. Figure [TJb). Let 




Figure 3: (Color online) (a) An arbitrary triangle ABC in its initial undefomied state. Its three edges have 
length di, d2, and ds (b) The triangle ABC is deformed into A'B'C This stretching deformation can be 
captured either by the changes in the edge lengths: Adi, Ad2, and Ads, or by the strain tensor, which is 
assumed constant across the triangle. 




Figure 4: (Color online) An arbitrary triangle ABC and its three adjoining triangles in both the initial 
undeformed stated and the deformed state in a local coordinate system. This bending deformation can be 
captured either by the three dihedral angles 6i, 62, and 6*3 between ABC and its three neighbors or by the 
curvature tensor, which is assumed constant across the four triangles. 

the point density of the initial uniform lattice be pinit- For our purposes we may choose a mapping 
with radial symmetry. For a point (r, 6) in the uniform lattice, we transform r so its new position 
in the variable lattice is (f(r), 9). Then the local point density in the variable lattice is pmit/ {%^)- 
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In practice, we choose r to be a function form of r = r + ^ gi{r, ri, Wi, Si), where 

i 

gi(r,ri,Wi, Si) = — [arctanf -) — arctanf -)], 

Wi Wi Wi 

dgjjr, ri,Wi,Si) ^ Si 

dr (r — Tj)^ + wf ' 

and the index i labels the regions where local point density needs adjustment. In Figure [IJb) there 
are four such regions: the center region, the rim region, the region between the center and the rim, 
and the outer region. In the center region and the rim region the point densities are increased while 
in the two other larger regions the point densities are reduced. The overall average point density 
remains close to that of the uniform lattice before the density adjustment. As shown in Figure 
[51 the graph of ^ is a simple U-shaped curve centered around r^, and we can control its width 
through Wi and depth through Sj. 
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Figure 5: (Color online) ^ vs. r when rj = 50, Wi = 20, and Si = —100. 

The resulting nonuniform lattice defines the initial undeformed state of Figure |3] for each tri- 
angle in the nonuniform lattice. Thus the lattice positions defined by this map constitute the state 
of zero stretching energy. By construction, these positions lie in a plane, so that each triangle also 
has zero curvature energy as defined by Figure SI 

Creating the nonuniform lattice for a c-cone requires two additional steps. To simulate a c-cone 
with opening angle equal to 29o, we make a cut along the radial line 6 = 0, and map every lattice 
point through the transformation r' = f and 9' = sin{9o)9. The last step is to join the two free 
boundary lines. That is, we identify each point on the free radial boundary line with its counterpart 
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on the other free boundary line, so that the two lattice positions are constrained to occupy the same 
spatial position in the simulation. 

The constraining container rim and the pushing force are simulated in almost the same way 



as previously used in Ref. [1241 13611 . The rim is in the x — y plane and its shape is determined 



by the equation x"^ + y"^ = R^, where R is the radius of the container. We introduce an external 



potential to implement the geometric rim constraint ( a/x^ + — R^ + z'^ ^ for all points in the 
sheet. To implement such a constraint for a discrete lattice we must assure that every lattice point 
remains a distance of order a from the rim line. For numerical tractability we assure this by adding 
an external potential felt by all lattice points that maintain the required separation while having 
negligible effect on more distant points. We find empirically that a short range potential is 
adequate. More specifically, we implement the repulsive normal force from the rim by introducing 
a potential of the form 



where Cp is a constant, {xj,yj, zj) is the coordinate of the jth lattice point, and the summation is 
over the whole lattice. The value of Cp is chosen so that the shortest distance between the lattice 
points and the rim is approximately the local lattice spacing in the radial direction. The potential 
due to the central pushing force is 

f/ibrce(a;i, yi, zi) = -{zi + a)G{xi,yi)F . 

Here (xi, yi, zi) is the coordinate of the lattice point in the center, F is effectively the magnitude 
of the pushing force, and the function G{xi, yi) is given by 

G{x,,y{) = [{l + {x^/if){l + {yl/m-\ 

where is a constant of order 0.1a. This G{xi, yi) is introduced to make sure that when the sheet 
is being pushed the lattice point in the center does not stray away from the axis of the cylindrical 
container, i.e., (xi, yi) ~ (0, 0). 

The total energy of the system is the sum of the total elastic energy in the sheet and the potential 



energies due to the rim and the applied pushing force. The conjugate gradient algorithm [|46ll is used 



to minimize the total e nergy a s a function of the coordinates of all lattice points to get the final 



shape of the sheet 124, 



4711 . To verify that the energy of the final configuration is at a global 



minimum, we can move each lattice point away from its equilibrium position by a random amount 
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in a random direction and see if the energy minimization process will bring the system back to its 
original state; the magnitudes of the random displacements introduced are usually much less than 
the local lattice spacing. We can also check how sensitive the final configuration is to the starting 
configuration. A high sensitivity generally signals the thickness of the sheet is too small for the 
current lattice to simulate and the results are likely unreliable. 

This model faithfully represents continuum sheets provided that the radii of curvature are every- 
where much larger than the local lattice spacing. This limitation restricts the values of deflection e 
to be less than or equal to 0.15 in our simulations. Ref. ll44n reports further details about this simu- 
lation technique. The simulation program using the nonuniform lattice has been validated to show 
uniform elastic behavior for states of planar stress and for cylinders. It has been validate d ag ainst 
uniform lattices for d-cones of thickness h that can be simulated by both methods. Ref. ll44ll also 
provides further simulated sheets, as well as our relaxation protocol and timing information. 

Our computer program is sequential and the bottleneck of the simulation is the CPU speed. A 
typical d-cone or c-cone is simulated with 67, 951 lattice points, and it usually takes more than 
one month to finish the energy minimization process on a 2.67 GHz Intel Core i7 processor. Each 
processor usually runs only one instance of the program at a time, but we have up to 40 processors 
to run multiple instances with different parameters simultaneously. 

Once the final configuration of a d-cone is obtained, the curvature tensor in each triangle can 
be determined using the procedure stated earlier. We then solve the characteristic equation of 
each curvature tensor to get the principal curvatures and Cee- To find the radial profile of 
Crr, we pick a circular sector with angle 0.05 and for each triangle within this sector, its Cj-r is 
associated with the radial coordinate of its center. As shown later, Crr reaches its maximum at 
the rim where the interaction between the sheet and the container is the strongest. Within each 
sector we average the four largest values of Crr near the peak as the rim value of Crr- The four 
corresponding Cqo are averaged to get the rim value of Cge. The rim curvatures depend on the 
angular separation between the center line of the circular sector and that of the buckled region. 
The angular separation is 57r/6 for the rim profiles of Crr and rim values of \Crr/Cgg \ reported in 
subsequent sections. The estimated percentage uncertainty of the reported rim values of \Crr/Cgg\ 
due to this angular dependence is about 10% for the whole range of thickness considered. A better 
approach is to average the rim Crr and Cgg across multiple sectors, but this will have no noticeable 
impact on our results. Our tests also show that different initial configurations will cause \Crr/Cgg\ 
at the rim to change by less than 4% for the range of thickness considered here. If we assume these 
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two sources of uncertainties are independent, the total percentage uncertainty for the reported rim 
values of ICrr/Ce^l should be about 11%. The results for c-cones are determined in the same way 
and have the same level of uncertainties as in the d-cone data. 

ni. NUMERICAL RESULTS 

Figured shows \Crr/Cee\ at the rim of d-cones versus relative thickness h/R for two different 
values of e. Unless explicitly stated otherwise, R and Young's modulus Y are assumed to be 
constant. For both e = 0.10 and e = 0.15, \C„-/Ceg\ goes below 1 as h/R is sufficiently thin, 
and it does not show any sign of leveling off as it reaches as low as 0.76 in the thinnest sheets 
simulated. More strikingly, for each fixed e, \Crr/Cee\ scales as [h/ RY^'^. This clearly contradicts 
the previous observation of vanishing mean curvature which requires \Crr/CeQ\ to stay at 1 when 
the sheet gets very thin. To resolve this contradiction, we need a better understanding of how 
\Crr/Cee\ at the rim responds to changes in h and F. However, in d-cones there is a one-to- 
one correspondence between h and F for a fixed e. To gain more flexibility and to explore the 
generality of this feature, we study d-cone's close relative, the c-cone. 




Figure 6: (Color online) The ratio of the two principal curvatures at the rim of a d-cone as a function of the 
relative thickness of the sheet for e = 0.10 and e = 0.15. It clearly shows that the ratio is less than one for 
very thin sheets and the ratio keeps decreasing as the sheet gets thinner. The slopes of the fitted lines are 
0.36 and 0.34 for e = 0.10 and e = 0.15 respectively. 

In c-cones, for a fixed e = 0.10, while h and F are changed independently we find that 
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\C„/C0e\ at the rim scales as {R/hY^'^F / {YB?), as shown in Figure |7l This scaling law is justi- 
fied in Section IV. In d-cones, even though h and F are interdependent, the same scaling law also 
holds. Moreover, Figure |8] demonstrates that for the same set of h, F, and e the radial profile of 
normalized C,.,. near the rim is the same in the d-cone as in the c-cone. These two findings together 
suggest very strongly that the behavior of Crr in the rim region of a d-cone is identical to that in a 
c-cone, and should be explained by the same mechanism. 



O 
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Figure 7: (Color online) \Crr/Ceg\ at the rim as a function of {R/hf/'^F/{YR^) for both c-cones and 
d-cones. h and F we. changed independently in the c-cone data. For a specific h, if we denote F^ as the 
center force needed in the d-cone to give the e of the c-cone, namely 0.1, the pushing forces used in the 
corresponding c-cone simulations may vary from 0.1 2Fd to l.SF^. 



To shed some light on the underlying mechanism, we investigate the detailed deformation of the 
sheets near the supporting rim, as sketched in Figure Eb). Let us denote the maximum deviation 
from the straight radial line as d and the width of the deformation as h. Figure |9] shows that d 
has a linear response to F and scales as {R/h)^^'^F/ (YR). From Figure [TOl we can see that b is 
independent of F and scales asVhR. It's worth noting that for the same h and within numerical 
accuracy b is exactly the same in a d-cone as in a c-cone. It scales in the same way as the width 
of a Pogorelov rin g rid ge formed by pushing a convex thin shell with a large inward concentrated 



normal force [|32 



371]. Our scaling arguments in the next section will closely follow how the 



scaling properties of the Pogorelov ring ridge are derived. 
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Figure 8: (Color online) Radial profiles of normalized in both a c-cone and a d-cone. The c-cone and 
the d-cone have the same e, h, R, and F. Each profile is nomialized by dividing the vertical scale by its 
corresponding peak value. There is no normalization for the horizontal scale. As stated in Section II, for the 
d-cone the angular separation between the buckled region and the radial sector used to generate the radial 
profile is 57r/6. 

X 10"'' 




Figure 9: (Color online) The maximum deviation d of the local deformation near the rim as a function of 
{R/hY^'^F / (YR). The constant container radius R is set to be the unit length. In the c-cone data h and F 
are changed independently as in Figure |7] 
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Figure 10: (Color online) The width b of the local deformation near the rim as a function of VhR. More 
specifically for the data shown here, width b is the fuU-width-at-half-maximum of the radial profile of Crr 
as shown in Figure [8] The constant container radius R is set to be the unit length. Again, h and F are 
changed independently in the c-cone data as in Figure |7] This plot shows that b does not depend on F in 
c-cones. 



IV, SCALING ARGUMENTS AND ANALYTICAL SOLUTIONS 



1. Scaling arguments 

Scaling arguments can be constructed to determine the dependence of ratio \Crr/Cee\, the max- 
imum deviation d, and width b on h, R, and F. Near the rim of a c-cone or a d-cone, the energy 
contribution of the local deformation includes both bending and stretching energy. Let us first find 
how each of them scales with d and b. 

To evaluate the bending energy component, we start with the principal curvatures. We denote 
quantities unperturbed by the rim force by an overbar, e.g., Cee. We denote changes induced by 
the the rim force by a A, e.g., ACgg. Let (p be the deviation of the tangent vector away from the 
original direction along a meridian as shown in Figure |2tb). In the deformed region ~ and 
~ (p/b ~ d/b^ . We also have Cgg ~ l/R, and radius of the azimuthal curvature 
Rgg = I /Ceo ~ R- Due to the local deformation, ^Ree ~ —d, so ACee ~ d/R^Q, or d/R^. Since 
b <^ Ras argued below, ACeg is much less than AC^r and can be safely ignored. Locally there is a 
contribution to the kg part of the bending energy of Equation |2] The change in this energy induced 
by the rim force is of order KcACrrCeo. However, this kq energy must integrate to zero because 
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of the Gauss-Bonnet theorem. Thus the bending energy UB,rim due to the local deformation near 
the rim (over an area ~ hR) is: 

UB.rim ~ nCrr'^hR ~ nRd^ /h\ (3) 



For the stretching component, d/R and 'jrr is negligible here pill . Thus the stretching 

energy contribution near the rim is: 

Us,r^m ~ hY^ee^bR ~ hYbd^R. (4) 

Minimizing U B,rim + Us,rim gives b ~ \/hR. Plugging this back into Equations [3] and H the 
total elastic energy is 



Us^rim + Us,rim ~ hYd^^Jh/R. (5) 

Taking the derivative of the total elastic energy with respect to d and equating the result to the 
pushing force F, we find d ~ {R/hf/^F/{YR). Thus Crr ~ d/b"^ ~ {R/hf/^F/{YR^), and 
\Crr/Cee\ ~ (R/hY^'^F/ (YR^). All the scaling relations obtained here agree with the numerical 
results presented in the previous section. 



2. Analytical solutions for c-cones 



The deformation of shells of revolution under symmetrical loading is a classical problem 141 
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490. The governing Foppl-von Karman equations [1320 are 4th order nonlinear differential 
equations. The analysis is greatly simplified by treating the limiting regime of weak loading, so 
that the equations can be linearized in the deformation. From the linearized solution, we may 
then verify that for loading forces of interest, the linearized treatment is completely valid for 
asymptotically thin sheets. 

For a c-cone, assuming that the supporting container's rim is infinitely hard (i.e. the range of 
the normal force from the rim is close to zero), satisfies the following differential equations as 
shown in Appendix B: 

.rfV d^(P ^ 1-^ forr <i?csc(eo) 



dr^ dr^ 



for r > _Rcsc( 



'0 



where is half of the underlying cone's opening angle, Ai = 12(1 — z/^) cot'^ (Oq) / h"^ , and A2 
Fsec\eo)/{27rYh). 
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These equations have closed form solutions in terms of Bessel functions. For r > i?csc 



the equation is homogeneous and its general solution is [|41 
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--B,[bei{0 + ^ber'iO] + B2[ber{C) - ^bei'iQ] 



(7) 



where C = 2^12(1 — 1^2^7.2 tan^(6'o) / h"^, and a prime is differentiation with respect to The ber, 
bei, ker, and kei functions are known as the Thomson or Kelvin functions [iJ 5^. Evidently 
goes to zero for large r. However, bei and ber diverges there, and are linearly independent. Thus 
their coefficients Bi and B2 must vanish. The solution for the r < R csc(6'o) has an extra term: 



-.B,[betiC) + ^ber'iC)] + B,[beriC) - ^bet'iQ] 
+ B,[kei{Q + -^ker'iO] + B,[ker{C) - ^A;ez'(C)] 



(8) 



Even though kei{C) and ker{Q are divergent at = 0, the kei and ker terms are needed to 
balance the A2/r term that also diverges at the apex. However, kei{() and ker{() decreases almost 
exponentially as increases, so these two terms are negligible near the rim. For our purpose 
of getting the radial profile of Crr near the rim, the term —A2/r can also be ignored because 
it changes in length scale ~ R, which is much larger than the width of the rim region. The 
contribution of this term to Cj-r at the rim is vanishingly small compared with the €„■ we got from 
numerical simulations and scaling arguments: 



d 
dr 



~A/r)/C^ 



A, 



YR^ 



. h ^ 
R' 



3/2 



(9) 



We conclude from the preceding reasoning that it is sufficient to specify S3, i?4, B^, and 
Bq for the inner and outer regions. To determine these four coefficients, we may use the four 
matching conditions applicable at the forcing point r = Rcsc(9q): is continuous and equal 
to zero , the curvature dcjy/dr is also continuous, and there is a jump for d^cp/d^r which equals 
— F/(27r sin(6'o)fi:-R). The last condition is due to the assumption that the rim of the supporting 
container is infinitely hard. Any localized force on an elastic sheet produces a discontinuity of 
curvature derivative of this type [|32||. 

Finally, we can compare the radial profile of C^r from analytical solution with that from nu- 
merical simulation. The overall excellent agreement between them, as shown in Figure [TTl gives 
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convincing evidence that our numerical results are valid. It should be stressed that there is no nor- 
malization in either axis. The 15% difference between the peak values are primarily due to three 
factors. First, as noted in Section II the percentage uncertainty of the rim curvatures from simu- 
lation is about 11%. Second, the analytical solution assumes an infinitely sharp container edge, 
but the force range of the normal force used in the simulation is close to 7% of the fuU-width- 
at-half-maximum (FWHM) of the peak for this specific c-cone. Third, the local lattice spacing 
in the radial direction is also about 7% of the FWHM. The last two factors cause the simulated 
curve to have a rounded peak. The effect of the finite thickness of the sheet is negligible here since 
the thickness is less than a tenth of the normal force range. Finite size of the simulation can also 
influence the curvature profile. 

In addition to the stronger peak, the analytical solution shown in Figure [TT]has a 10% stronger 
dip on either side. We believe this is due to the local compensation of Gaussian curvature [13 611 . 
which requires the integral under the curve to be zero. So if the analytical solution has extra area 
under the peak, it must have extra negative area at the dips. 

Over the full range of the c-cone sheet thickness covered in our simulation, the peak values 
of the Crr at the rim from simulation is lower than that given by the analytical solutions by be- 
tween 10% and 15%. We should expect a similar level of discrepancy between the simulation and 
analytical solutions for the d-cone. Thus the lowest \Crr/Cee\ value achieved for d-cones in our 
simulations may increase from 0.76 to a value as high as 0.87, which is much closer to 1, but this 
level of discrepancy should have no material impact on the scaling relationship between | Cee \ , 
h, and F, which is the much stronger evidence that \Crr/Cee\ at the rim will drop below one and 
keep decreasing as the thickness of the sheet approaches zero. 



V. DISCUSSION AND CONCLUSION 



In this paper we have shown numerically that contrary to previous claims, the mean curvature 
at the rim in a d-cone does not vanish as the thickness of the sheet goes to zero. This vanishing 
requires that | Cj-r/ Cge \ goes to 1 . However, in the range we studied, | €„/ Cee I at the rim appears to 
vary as {h/RY^^. More generally, in both c-cones and d-cones, \Crr/Cee\ ~ {R/h)^/'^F/{YE?). 
These identical scaling laws and the similarity of the radial profiles of €„ in d-cones and c-cones 
suggest that the core region of a d-cone has no influence on how the rim region reacts to the normal 
rim force pressure. The core region only affects the amplitude of the rim force pressure. 
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Figure 1 1 : (Color online) Radial profiles of Crr in a c-cone from both numerical simulations and analytical 
solutions. Crr is in units of l/R. The peak value from simulations is about 15% lower than that from 
analytical solutions. The sheet here is the thinnest used in simulations with relative thickness h/R = 
0.00029. 



This paper does not attempt to determine the right scaling law for the pushing force F in a 
d-cone. Simply combining the derived h~^^'^F scaling law with the numerically observed h^^^ 
scaling law for \Crr/Cge\ will give F ~ Yh^ / R{h/ R)~^/^ . There is another proposed functional 
form for F: F Y h? / R\n.{RJ R^) , where Rp is radius of the sheet, and Rc ~ h^/^R^/^ is the 
radius of the core region |l, 25]. The first force scaling is asymptotically much stronger than 
the second one with the logarithmic term. However, within the dynamic range covered in our 
simulations, the fit for these two functional forms are equally good. It should also be mentioned 
that some researchers have expressed doubts on the arguments leading to the second functional 



form[|24, 



2711 ■ Our work in progress aims to resolve this issue. 
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Appendices 



Appendix A: Formulas for strain and curvature tensors 
1. Formulas for the strain tensor 

This subsection derives the expression for strain tensor 7 in terms of the changes of edge lengths 
Adi, Ad2 and Ad^ in an arbitrary triangle ABC as shown in Figure Ha) and|3lb). 

In FigureOa), let the coordinates of A, B, and C be (0, 0), {xb, Ub), and (xc, 0), respectively. 
When triangle ABC is under strain, as shown in Figure Hb), denote the changes of x^, Ub, and xc 
as Axb, Ays, and Axc, and the new coordinates of A, B, and C are (0, 0), {xb+Axb, yB+AyB), 
and {xc + AxcO). For a general point in the triangle, let us denote its displacement 

as {ux,Uy), so its coordinate due to deformation is (x + Ux,y + Uy). Under the assumption of 
infinitesimal constant strain across the triangle, and ^ should be constants. We 

^ ox ^ oy ^ ax ^ oy 

can also expand the displacement of points B and C in terms of these partial derivatives: 

Axb = -tt-xb + -tt-Ub + higher order terms , (Al) 

ox oy 

du du 

Ays = -tt^xb + -TT^yB + higher order terms , (A2) 
dx dy 

du 

Axc = -TT^Xc + higher order terms , (A3) 

dx 

du 

Aye = -TT^^c + higher order terms = . (A4) 

dx 

Ignoring the higher order terms, we can solve Equations lAlJ LA 2 1 A3 [ and I A4I for and 

dUy _ 

dy ■ 

du^ Axc 



(A5) 



dx Xc 

du^ Axb xbAxc 

= , (A6) 

dy Vb VbXc 



dUy 

dx 

duy AyB 



dy Vb 



0, (A7) 

(A8) 
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Then we can easily get the strain tensor elements [l32ll : 

du^ Axc 



^yy 



dx 
dun, 



dy 

Ixy = 7ya; 



Vb 

1 ,dux , duy 1 Axb xbAxc. 



2'dy 



+ 



dx 



2' Vb 



VbXc 



(A9) 
(AlO) 
(All) 



or in matrix form: 







Axb 


lyy 


= S 


^Vb 


Ixy 




Axc 



(A12) 



where 



S 









1 

xc 





1 

ys 





1 





XB 


2yB 


"iyBxc 



(A13) 



We have determined the elements of 7 in terms of Axb, AyB, and Axc, but it is much more 
efficient to compute Adi, Ad2, and Ad, during simulating, so let us find the expression for Axb, 
AyB, and Axc in terms of Adi, Ad2, and Ad,, and then express 7 in terms of Adi, Ad2, and 
Ad,. 

The change of length for edge Ai? is: 



Arfi = ^/{^^^^TA;^^^yV(y^TA^ -\/xl + y% 
A , Vb 



xl + vl 



Axb + 



--AyB + higher order terms 



Similarly, we can express Ad2 and Ad, in terms of Axb, AyB, and Axc'- 

Xb- Xc , ^ _^ A „ X , 



(A14) 



Ad2 = - 

V{xb - xcf + y\ 
+ higher order terms 

Ad, = Axc . 



{Axb - Axc) 



^ {xb - xcY + y\ 



Ay_ 



B 



(A15) 
(A16) 



Let us ignore the higher order terms and rewrite Equations lA 1 41 IaTSI and IA 161 in matrix form: 



Adi 




Axb 


Ad2 


= G 


^VB 


Ad, 




Axc 



(A17) 
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where 



G 



XB-XC 



Vb 
xg-XB 







Vb 



y/ {xB-xc)'^+y% \J{xB-xc?-+y'B \J {xB-xc)^+y% 

1 



Combining Equations IA12I and IA17I 



'fxx 






lyy 






Ixy 




Ads 



(A18) 



(A19) 



where G^^ means the inverse of G. Since both S and G depend only on xb, Vb, and xc, the 
matrix SG^^ is determined by the initial geometry of the triangle and can be calculated at program 
initialization. 



2. Formulas for the curvature tensor 

To calculate the curvature tensor on an arbitrary triangle, e.g. triangle ABG in Figure ID we 
can fit the coordinates of the six vertices of its three adjoining triangles to the following function 
[0,3: 

Zi = ai + a2Xi + a^Vi + a^xl + a^Xiyi + agi/f , i = A,... ,F (A20) 

where (xj, y^, Zi) are coordinates of the vertices in a local coordinate system. In this system, the z 
axis is perpendicular to ABG and its origin is at the center of ABG. These choices ensure that 02 



and 03 are negligible. Then the curvature tensor elements are determined through [[19 



m 



Gxx — 2a4 , Gxy — 05 ) — • (-^-21) 

If we ignore the changes in the Xi and yi, the coefficient matrix in Equation IA20I will stay the 
same during the simulation, which means for each triangle we only need to invert the coefficient 
matrix once at program initialization. Under the assumption that both bending and stretching are 
infinitesimal, this simplification will result in a second order error in the curvature tensor. 

In this local coordinate system, za, zb, and zc are all zero by the choice of the z axis. The 
other three nonzero z coordinates zd, ze, and zp can be determined from the three dihedral angles 
between ABG and its three adjoining neighbors. For example, zd ~ dD,AcOi, where dD,Ac is the 
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distance from point D to edge AC. So, alternatively, we can also determine the curvature tensor 
from the three dihedral angles 9i, 02, and 6*3. 

This method works for both uniform and variable lattices. 



Appendix B: Analytical solutions for c-cones 




Figure 12: (Color online) Edge view of an element along the meridian in a c-cone. W and T are the axial 
and radial forces per unit length, p is the radial distance from the axis, and r is the meridional length 
(equivalent to the radial distance from the apex in the material coordinate). For an ideal c-cone, the normal 
pressure p is zero everywhere except at the rim. F and Oq are the central pushing force and the half opening 
angle, respectively. 



Ref . ll42h provides detailed information on the general theory of symmetrically loaded shells of 
revolution, including conical shells. Its derivation assumes that the deformation is small relative 
to the size of the structure, but may be comparable to the thickness. In this regime, both bending 
and in-plane stretching may occur and thus need to be considered together. However, for many 
problems only the force and moment balance of the undistorted element is needed, and the result- 
ing equations are normally linear. Here, we will simply quote the relevant equations presented 
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in Ref. [14211 and use them to find the specific equations for c-cones assuming that the supporting 
container's rim is infinitely hard. 

As shown in Figure \T2l p, W, and T are the pressure, axial and radial forces per unit length, 
respectively, p is the radial distance from the axis and it is related to r through p = r sin(6'o). In a 



conical shell these variable and satisfy the following equations according to Ref. [|42l1 : 

Wp = - J ppdp + V , (Bl) 

«:sin(^o)[p^{-^(p0)}]+Tpcot(^o) = -y Ppdp + V, (B2) 

sm{eo)[p^{-^{p'T)}] - hY^cot{9o) = cos{eo)[^ {pp') - - / ppdp + -] , (B3) 
dp pdp dp p J p 

where V is a constant of integration. It is worth noting that these equations are derived for sheets 
whose unstressed (or undeformed) state has curvature, but we find that they are also valid when 
the unstressed state is flat. 

Since we assume the supporting container's rim is infinitely hard, using the balance of force it 
is straightforward to determine that in c-cones W has the following functional form: 

W = ^H{R-p) (B4) 
Znp 

where H{x) is the Heaviside step function. Combining Equations IB 1 1 and lB4l we can get: 

P=7^S{p-R),mdV=^, (B5) 

where S{x) is the Dirac delta function. Plugging this into Equations IB2I and IB3I and replacing p 
with r sin (^o) yield: 

d'(j) d(j) (j) ^ Trcos(^o) forr <i?csc(^o) 

+ -J \ = \ , i^oj 



dr'^ dr r 



for r > _Rcsc( 



^0) 



^ d\Tr) ^ d{Tr) (Tr) hY<P cos{9o) _ J/.-ffgi) for r < i?csc(^o) ^^^^ 
' ' ^i^'(^o) [o forr>i?csc(^o) ' 

Equation IB6I can be rewritten to get an explicit expression for Tr: 

d'^d dd 6 \^^^ forr < i?csc(^o) 



dr"^ dr r ' 
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for r > i?csc( 



'0) 



Plugging Equation IB8I into Equation IB71 we can get a forth order partial differential equation 
for 0: 



which is equivalent to Equation [6l 
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